TLS

Preamble

Dependencies

Code

library(scran)
library(scater)
library(scuttle)
library(ggplot2)
library(pheatmap)
library(patchwork)

Loading

Code

spe <- readRDS("../outs/01-spe.rds")

Setup

Wrangling

Code

# add tissue types
spe$TissueType <- ifelse(
    spe$TumorType == "ccRCC", 
    yes = "kid", no = "lun")
# exclude unassigned & 
# spots tagged for exclusion
nan <- is.na(spe$anno1) | is.na(spe$anno2)
rmv <- grepl("EXCL", spe$anno1)
sub <- spe[, !(nan | rmv)]
# subset regions of interest
ids <- c("INFL", "TLS", "LN")
sub <- sub[, sub$anno1 %in% ids]
sub$anno1 <- droplevels(sub$anno1)
# subset samples of interest
ids <- c("B04_17776", "B06_24137", "B06_24784")
sub <- sub[, sub$sample_id %in% ids]
sub$sample_id <- droplevels(sub$sample_id)
# simplify annotations
sub$anno3 <- as.character(sub$anno2)
sub$anno3[grep(".*ETLS", sub$anno3)] <- "E_TLS"
sub$anno3[grep(".*MTLS", sub$anno3)] <- "M_TLS"
sub$anno3 <- factor(sub$anno3, exclude = NULL)
table(sub$sample_id, sub$anno3)

           
            E_TLS INFL  LN M_TLS
  B04_17776    68  140   0    94
  B06_24137   235   62  33   166
  B06_24784    89  173   0    53

Code

table(sub$anno1, sub$anno3)

      
       E_TLS INFL  LN M_TLS
  INFL     1  375   0     5
  TLS    391    0   0   308
  LN       0    0  33     0

Filtering

Code

# keep features detected in at least 20 spots in any sample
# and spots with at least 200 detected features overall
idx <- split(seq(ncol(sub)), sub$sample_id)
gs <- sapply(idx, \(.) {
    y <- counts(sub[, .]) > 0
    rowSums(y) >= 20
})
fil <- sub[rowAnys(gs), ]
fil <- fil[, colSums(counts(fil) > 0) >= 200]
cbind(spe = dim(spe), sub = dim(sub), fil = dim(fil))

       spe   sub  fil
[1,] 16643 16643 8814
[2,] 27595  1113 1065

Analysis

Code

# split by sample
idx <- split(seq(ncol(sub)), sub$sample_id)
lys <- lapply(idx, \(.) sub[, .])
# feature selection & PCA
lys <- lapply(lys, \(.) {
  tbl <- modelGeneVar(.) 
  hvg <- rowData(.)$hvg <- tbl$bio > 0
  runPCA(., subset_row = hvg) 
})
# within TLS only
tls <- lapply(lys, \(.) {
  . <- .[, .$anno1 == "TLS"]
  .$anno3 <- droplevels(.$anno3)
  tbl <- modelGeneVar(.) 
  hvg <- rowData(.)$hvg <- tbl$bio > 0
  runPCA(., subset_row = hvg) 
})

Code

thm <- list(
  coord_equal(),
  theme(legend.key.size = unit(0.5, "lines")),
  guides(col = guide_legend(override.aes = list(alpha = 1, size = 2))))
pal <- hcl.colors(nlevels(sub$anno3), "Set 2")

lapply(lys, plotPCA, color_by = "anno3") |>
  wrap_plots() + plot_layout(guides = "collect") & 
  thm & scale_color_manual(values = pal, drop = FALSE)

Code

lapply(tls, plotPCA, color_by = "anno3") |>
  wrap_plots() + plot_layout(guides = "collect") & 
  thm & scale_color_manual(values = c("red", "blue"))

Code

var <- sapply(tls, \(.) {
  var <- modelGeneVar(.)$bio
  setNames(var, rownames(.))
})
# top variables across all sample
o <- order(rowMeans(var), decreasing = TRUE)
(sel <- rownames(var)[head(o, 12)])

 [1] "IGHG1"   "IGKV4.1" "IGHG2"   "IGLC1"   "IGLV3.1" "IGHG3"   "IGHA1"  
 [8] "IGHG4"   "CLU"     "IGHJ6"   "XBP1"    "H4C4"   

Code

lapply(tls, \(sce) {
  lapply(sel, \(g) plotPCA(sce, color_by = g)) |>
    wrap_plots() + plot_layout(nrow = 3) & coord_equal()
}) 

$B04_17776

$B06_24137

$B06_24784

Code

# top variables by sample
sel <- apply(var, 2, \(.) {
  o <- order(., decreasing = TRUE)
  names(.)[head(o, 20)]
}, simplify = FALSE)
lapply(names(tls), \(.) {
  x <- tls[[.]][sel[[.]], ]
  cd <- data.frame(colData(x))
  cd <- cd[c("total", "anno3")]
  hm <- pheatmap(
    main = .,
    logcounts(x),
    scale = "none",
    annotation_col = cd,
    show_colnames = FALSE)
  print(hm)
})

[[1]]
NULL

[[2]]
NULL

[[3]]
NULL

Appendix

Code

sessionInfo()

R version 4.3.2 (2023-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.6 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=de_CH.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=de_CH.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=de_CH.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=de_CH.UTF-8 LC_IDENTIFICATION=C       

time zone: Europe/Zurich
tzcode source: system (glibc)

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] SpatialExperiment_1.10.0    patchwork_1.1.3            
 [3] pheatmap_1.0.12             scater_1.28.0              
 [5] ggplot2_3.4.3               scran_1.28.2               
 [7] scuttle_1.10.2              SingleCellExperiment_1.22.0
 [9] SummarizedExperiment_1.30.2 Biobase_2.60.0             
[11] GenomicRanges_1.52.0        GenomeInfoDb_1.36.3        
[13] IRanges_2.34.1              S4Vectors_0.38.1           
[15] BiocGenerics_0.46.0         MatrixGenerics_1.12.3      
[17] matrixStats_1.0.0          

loaded via a namespace (and not attached):
 [1] bitops_1.0-7              gridExtra_2.3            
 [3] rlang_1.1.1               magrittr_2.0.3           
 [5] compiler_4.3.2            DelayedMatrixStats_1.22.6
 [7] vctrs_0.6.3               pkgconfig_2.0.3          
 [9] crayon_1.5.2              fastmap_1.1.1            
[11] magick_2.7.5              XVector_0.40.0           
[13] labeling_0.4.3            utf8_1.2.3               
[15] rmarkdown_2.24            ggbeeswarm_0.7.2         
[17] xfun_0.40                 bluster_1.10.0           
[19] zlibbioc_1.46.0           beachmat_2.16.0          
[21] jsonlite_1.8.7            rhdf5filters_1.12.1      
[23] DelayedArray_0.26.7       Rhdf5lib_1.22.1          
[25] BiocParallel_1.34.2       irlba_2.3.5.1            
[27] parallel_4.3.2            cluster_2.1.5            
[29] R6_2.5.1                  RColorBrewer_1.1-3       
[31] limma_3.56.2              Rcpp_1.0.11              
[33] knitr_1.44                R.utils_2.12.2           
[35] Matrix_1.6-1              igraph_1.5.1             
[37] tidyselect_1.2.0          rstudioapi_0.15.0        
[39] abind_1.4-5               yaml_2.3.7               
[41] viridis_0.6.4             codetools_0.2-19         
[43] lattice_0.22-5            tibble_3.2.1             
[45] withr_2.5.1               evaluate_0.21            
[47] pillar_1.9.0              generics_0.1.3           
[49] RCurl_1.98-1.12           sparseMatrixStats_1.12.2 
[51] munsell_0.5.0             scales_1.2.1             
[53] glue_1.6.2                metapod_1.8.0            
[55] tools_4.3.2               BiocNeighbors_1.18.0     
[57] ScaledMatrix_1.8.1        locfit_1.5-9.8           
[59] cowplot_1.1.1             rhdf5_2.44.0             
[61] grid_4.3.2                DropletUtils_1.20.0      
[63] edgeR_3.42.4              colorspace_2.1-0         
[65] GenomeInfoDbData_1.2.10   beeswarm_0.4.0           
[67] BiocSingular_1.16.0       HDF5Array_1.28.1         
[69] vipor_0.4.5               cli_3.6.1                
[71] rsvd_1.0.5                fansi_1.0.4              
[73] S4Arrays_1.0.6            viridisLite_0.4.2        
[75] dplyr_1.1.3               gtable_0.3.4             
[77] R.methodsS3_1.8.2         digest_0.6.33            
[79] ggrepel_0.9.3             dqrng_0.3.1              
[81] farver_2.1.1              rjson_0.2.21             
[83] htmlwidgets_1.6.2         htmltools_0.5.6          
[85] R.oo_1.25.0               lifecycle_1.0.3          
[87] statmod_1.5.0